      SUBROUTINE TOPCENXYZ(GLAT,GLONG,HEIGHT,RADEQTR,FLAT,XYZ)
      IMPLICIT REAL*8 (A-H,O-Z)
C
C FOR AN ELLIPSOID OF REVOLUTION, DO THE FOLLOWING --
C 
C    GIVEN: 1. A POINT OF INTEREST LYING A KNOWN DISTANCE ABOVE/BELOW A
C                 POINT ON THE ELLIPSOID. THE POINT ON THE ELLIPSOID HAS
C                 KNOWN GEODETIC LATITUDE AND LONGITUDE;
C           2. ELLIPSOID RADIUS AT THE EQUATOR;
C           3. FLATTENING COEFFICIENT OF THE ELLIPSOID;
C
C    FIND:  THE CARTESIAN COORDINATES OF THE POINT. 
C
C           THE AXES TO WHICH THESE COORDINATES ARE REFERENCED HAVE 
C           THEIR ORIGIN AT THE ELLIPSOID'S ORIGIN, THE XY PLANE 
C           COINCIDENT WITH THE ELLIPSOID EQUATOR, THE Z AXIS AS THE 
C           ELLIPSOID AXIS OF REVOLUTION(IE, THE AXIS ABOUT AN ELLIPSE 
C           WAS ROTATED TO PRODUCE THE ELLIPSOID). POSITIVE LATITUDE IS
C           MEASURED TOWARD +Z AND POSITIVE LONGITUDE IS MEASURED FROM 
C           THE XZ PLANE TOWARD THE YZ PLANE.
C
C
C VARIABLE DIM TYPE I/O DESCRIPTION
C -------- --- ---- --- -----------
C
C GLAT      1  R*8   I  THE GEODETIC LATITUDE OF A POINT ON THE 
C                       SURFACE BELOW THE POINT OF INTEREST. IN RADIANS.
C
C                       A LINE FROM THE POINT OF INTEREST TO THE 
C                       SURFACE IS NORMAL TO THE SURFACE AT LATITUDE 
C                       GLAT.
C
C                       GEODETIC LATITUDE IS THE ANGLE BETWEEN A LINE
C                       NORMAL TO THE LOCAL HORIZON AND THE EQUATOR.
C                       GEODETIC LATITUDE INVOLVES FLATTENING.
C
C GLONG     1  R*8   I  THE LONGITUDE OF THE POINT OF INTEREST.
C                       IN RADIANS.
C
C HEIGHT    1  R*8   I  THE HEIGHT OF THE POINT OF INTEREST ABOVE THE 
C                       ELLIPSOID. MEASURED NORMAL TO THE LOCAL HORIZON.
C                       ANY LENGTH UNIT MAY BE USED.
C
C RADEQTR   1  R*8   I  THE RADIUS AT THE EQUATOR OF THE ELLIPSOID.
C                       HEIGHT AND RADEQTR MUST HAVE THE SAME LENGTH 
C                       UNIT.
C
C                       THE EQUATOR LIES IN THE 'FAT' PART OF THE 
C                       ELLIPSOID, NORMAL TO THE AXIS OF REVOLUTION.
C
C                       IF RADEQTR IS ZERO OR NEGATIVE, EARTH RADIUS IN
C                       KILOMETERS IS USED.
C
C FLAT      1  R*8   I  THE FLATTENING COEFFICIENT. UNITLESS.
C
C                       SPHERE HAS FLAT=0.0
C
C                       IF FLAT < 0.0 IS ENTERED, OBLATE EARTH IS USED
C                       AS DEFINED BY CONST(59). CONST IS A FUNCTION
C                       SUBROUTINE.
C
C XYZ       3  R*8   O  THE CARTESIAN COORDINATES OF THE POINT.
C                       UNITS ARE THE SAME AS FOR RADEQTR AND HEIGHT.
C
C***********************************************************************
C
C BY C PETRUZZO, 10/83.
C
C***********************************************************************
C
C  EQUATIONS ARE FROM THE BOOK: METHODS OF ORBIT DETERMINATION, BY
C  P.R. ESCOBAL, 1983 REPRINT(WITH CORRECTIONS) OF 1965 EDITION,
C  PAGES 399 AND 400. 
C
C
      REAL*8 XYZ(3)
      INTEGER INIT/1/
      REAL*8 HALFPI/ 1.570796326794897D0 /
      REAL*8 DEGRAD/ 57.29577951308232D0 /
      LOGICAL SPHERE
C
      IBUG=0
      LUBUG=19
C
      IF(INIT.EQ.1) THEN
        INIT=0
        POLE = HALFPI - 1.D-3/DEGRAD
        REARTH = CONST(53)
        FLATCONST = CONST(59)
        END IF
C
      IF(IBUG.NE.0) WRITE(LUBUG,9001) GLAT*DEGRAD,GLONG*DEGRAD,HEIGHT,
     *     RADEQTR,FLAT
 9001 FORMAT(' TOPCENXYZ DEBUG.'/,
     *       '    GLAT=',G16.8,'  GLONG=',G16.8,'  HEIGHT=',G16.8/,
     *       '    RADEQTR=',G16.8,'  FLAT=',G16.8)
C
C  SET INTERNAL VARIABLE CONTAINING THE FLATTENING COEFFICIENT.
      IF(FLAT.LT.0.D0) THEN
        FL=FLATCONST
      ELSE
        FL=FLAT
        END IF
      SPHERE = FL.EQ.0.D0
C
C  SET THE INTERNAL VARIABLES CONTAINING LATITUDES.
C  GDLAT IS GEODETIC, GCLAT IS GEOCENTRIC.
      GDLAT=GLAT
      IF(SPHERE) THEN
        GCLAT=GLAT
      ELSE
        IF(DABS(GLAT).LT.POLE) THEN
          GCLAT=DATAN( (1.D0-FL)**2 * DTAN(GLAT) )
        ELSE
          GCLAT=GLAT
          END IF
        END IF
C
C  GET RC, THE RADIUS FROM THE ORIGIN OF THE ELLIPSOID TO THE 
C  SURFACE POINT BELOW THE POINT OF INTEREST.
      AE=RADEQTR
      IF(SPHERE) THEN
        RC=AE
      ELSE
        TEMP=(2.D0 - FL) * FL
        RC = AE * DSQRT( (1.D0-TEMP)/(1.D0-TEMP*( DCOS(GCLAT)**2) ))
        END IF
C
      IF(IBUG.NE.0) WRITE(LUBUG,9002) GDLAT*DEGRAD,GCLAT*DEGRAD,RC
 9002 FORMAT(' TOPCENXYZ DEBUG.'/,
     *       '    GDLAT=',G16.8,'  GCLAT=',G16.8,'  RC=',G16.8)
C
C  GET R, THE RADIUS FROM THE ORIGIN TO THE POINT OF INTEREST.
      HS=HEIGHT
      IF(SPHERE) THEN
        TEMP=1.D0
      ELSE
        TEMP=DCOS(GDLAT-GCLAT)
        END IF
      IF(HS.NE.0.D0) THEN
        R = DSQRT( RC**2 + HS**2 + 2.D0*RC*HS*TEMP )
      ELSE
        R = RC
        END IF
C
C  GET DELGCLAT, THE DIFFERENCE BETWEEN THE GEOCENTRIC LATITUDE OF THE
C  POINT OF INTEREST AND THE GEOCENTRIC LATITUDE OF THE SURFACE POINT
C  BELOW THE POINT OF INTEREST.
      IF(SPHERE) THEN
        DELGCLAT = 0.D0
      ELSE
        DELGCLAT = DASIN( HS/R * DSIN(GDLAT-GCLAT) )
        END IF
C
C  GET SIN/COSINES NEEDED BELOW.
      THETA = GLONG
      DELTA = GCLAT + DELGCLAT
      COSDELTA = DCOS(DELTA)
      SINDELTA = DSIN(DELTA)
      COSTHETA = DCOS(THETA)
      SINTHETA = DSIN(THETA)
C
      IF(IBUG.NE.0) WRITE(LUBUG,9003) THETA*DEGRAD,DELTA*DEGRAD,R
 9003 FORMAT(' TOPCENXYZ DEBUG.'/,
     *       '    THETA=',G16.8,'  DELTA=',G16.8,'  R=',G16.8)
C
C  COMPUTE CARTESIAN POSITIONS
      XYZ(1) = R * COSDELTA * COSTHETA
      XYZ(2) = R * COSDELTA * SINTHETA
      XYZ(3) = R * SINDELTA
C
      IF(IBUG.NE.0) WRITE(LUBUG,9004) XYZ
 9004 FORMAT(' TOPCENXYZ DEBUG. XYZ=',3G16.8)
C
      RETURN
      END
